---
title: "CJPS replication"
author: "Olivier Bergeron-Boutin"
date: "24/11/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
```

# Packages and data

```{r}
# Packages
library(tidyverse)
library(ggplot2)
library(estimatr)
library(modelsummary)
library(ggtext)
library(DescTools)
library(extrafont)
library(cregg)
library(lmtest)
library(margins)
library(ggpubr)
library(fastDummies)
library(vtable)

# Data
load("data/can_covid_clean.RData")

# Conjoint formulas
f1 <- selected ~ Gender + Age + `Political experience` + Party + 
  `Economic aid policy` + `Lockdown policy` + `Legislative checks`

f2 <- selected ~ Gender + Age + `Political experience` + Party + 
  `Economic aid policy` + `Lockdown policy` + `Judicial checks`
```

## Custom ggplot themes

```{r}
# taken from Andrew Heiss' website
library(ggtext)
theme_custom <- function(base_size = 16, legend.position = "none"){
  theme_bw(base_size = base_size,
           base_family = "Fira Sans") %+replace%
  theme(legend.position = legend.position,
        panel.grid.minor = element_blank(),
        plot.title = element_markdown(face = "bold", size = rel(1.7)),
        plot.subtitle = element_markdown(face = "plain", size = rel(1.3)),
        axis.title = element_text(face = "bold"),
        axis.title.x = element_text(margin = margin(t = 10), hjust = 0),
        axis.title.y = element_text(margin = margin(r = 10), hjust = 1, angle = 90))
}

theme_cjoint <- function(){
  theme_bw(base_size = 18,
           base_family = "Fira Sans") %+replace%
    theme(legend.position = "none",
          axis.text.y = element_text(hjust=0),
          axis.title.x = element_text(size = 14))
}
```

## Custom functions

```{r}
source("functions/prepare_conjoint.R")
source("functions/plot_conjoint.R")
source("functions/conjoint_scenario.R")
source("functions/confints.R")
```

# Study 1

## Regression models

```{r}
# final dataset for regressions
# will be useful to run other models that would otherwise use a different
# sample due to listwise deletion 
survey_reg <- survey %>% 
  filter(partyid != "Other") %>% 
  drop_na(ea_2item, anxiety_scale, birth_decade, educ_4cat, 
                deprivation_scale, authority_scale_alt, partyid)

# Model only with anxiety
m1_anxiety <- lm(
  ea_2item ~ anxiety_scale,
  data = survey_reg
)

# formula for the fully-specified model
reg_formula <- ea_2item ~ anxiety_scale + birth_decade + educ_4cat +
  deprivation_scale + authority_scale_alt + partyid

# Fully-specified model
m2_anxiety <- lm(
  formula = reg_formula,
  data = survey_reg
)

# Covariate names for regression table
cm_m2 <- c(
  "anxiety_scale" = "COVID-related anxiety",
  "birth_decade1950s" = "Born in the 1950s",
  "birth_decade1960s" = "Born in the 1960s",
  "birth_decade1970s" = "Born in the 1970s",
  "birth_decade1980s" = "Born in the 1980s",
  "birth_decade1990s" = "Born since 1990",
  "deprivation_scale" = "Deprivation",
  "authority_scale_alt" = "Authoritarianism",
  "educ_4catCompleted high school" = "Completed high school",
  "educ_4catSome postsecondary" = "Some postsecondary",
  "educ_4catCollege graduate" = "College graduate",
  "partyidConservative" = "Conservative partisan",
  "partyidNDP" = "NDP partisan",
  "partyidGreen" = "Green partisan",
  "partyidNon-partisan" = "Non-partisan",
  "(Intercept)" = "Constant"
)

# Regression table with output to word document
# Robust SEs set using the "vcov" argument
modelsummary(list(m1_anxiety, m2_anxiety),
             coef_map = cm_m2,
             out = "tables/table_a3.docx",
             stars = c("*" = .05, "**" = 0.01, "***" = 0.001),
             estimate = "{estimate} ({std.error}){stars}",
             statistic = NULL,
             vcov = "HC0")
```


## Figure 1 (predicted values)

```{r}
m2_pred <- prediction(
  m2_anxiety, 
  at = list(anxiety_scale = c(0, 0.165, 0.33, 0.5, 0.66, 0.83, 1))
) %>% summary()

ggplot(m2_pred, aes(x = `at(anxiety_scale)`, y = Prediction, 
                    ymin = lower, ymax = upper)) +
  geom_point() +
  geom_ribbon(alpha = 0.4) +
  geom_line(alpha = 0.6, group = 1) +
  theme_custom() +
  labs(x = "COVID-related anxiety",
       y = "Predicted support for executive aggrandizement") +
  scale_y_continuous(limits = c(0.1, 0.55))

ggsave("figures/study 1/fig_1.png", height = 7, width = 9)
```

# Study 2

## Lockdown congruence

```{r}
# dummy = 1 if respondent agrees with hypothetical premier on lockdownns
# 0 otherwise
# "lockdown until fewer deaths" is treated as pro-lockdown position
survey_reg <- survey_reg %>% 
  mutate(
    resp_pro_lockdown = case_when(
      covid_lockdowns == "Lockdown ended now" ~ 0,
      covid_lockdowns %in% c("Lockdown until fewer deaths",
                             "Lockdown until vaccine") ~ 1
    ),
    exp1_pro_lockdown = ifelse(
      str_detect(exp1_condition, "lockdown"), 1, 0
    ),
    exp1_lockdown_congruence = ifelse(
      resp_pro_lockdown == exp1_pro_lockdown, 1, 0
    )
  )
```

## Regression models

```{r}
m3_anxiety <- lm(
  formula = update(reg_formula, exp1_outcome ~ . + exp1_lockdown_congruence),
  data = survey_reg
)

m4_anxiety <- lm(
  formula = update(reg_formula, exp1_outcome ~ . + exp1_lockdown_congruence +
                     anxiety_scale:exp1_lockdown_congruence),
  data = survey_reg
)

# Covariate names for regression table
cm_m4 <- c(
  "anxiety_scale" = "COVID-related anxiety",
  "exp1_lockdown_congruence" = "Lockdown preferences congruence",
  "anxiety_scale:exp1_lockdown_congruence" = "COVID-related anxiety*lockdown preferences congruence",
  "birth_decade1950s" = "Born in the 1950s",
  "birth_decade1960s" = "Born in the 1960s",
  "birth_decade1970s" = "Born in the 1970s",
  "birth_decade1980s" = "Born in the 1980s",
  "birth_decade1990s" = "Born since 1990",
  "deprivation_scale" = "Deprivation",
  "authority_scale_alt" = "Authoritarianism",
  "educ_4catCompleted high school" = "Completed high school",
  "educ_4catSome postsecondary" = "Some postsecondary",
  "educ_4catCollege graduate" = "College graduate",
  "partyidConservative" = "Conservative partisan",
  "partyidNDP" = "NDP partisan",
  "partyidGreen" = "Green partisan",
  "partyidNon-partisan" = "Non-partisan",
  "(Intercept)" = "Constant"
)

# Regression table with output to word document
# Robust SEs set using the "vcov" argument
modelsummary(list(m3_anxiety, m4_anxiety),
             coef_map = cm_m4,
             out = "tables/table_a4.docx",
             stars = c("*" = .05, "**" = 0.01, "***" = 0.001),
             estimate = "{estimate} ({std.error}){stars}",
             statistic = NULL,
             vcov = "HC0")
```

## Figure 2 (predicted values)

```{r}
m3_pred <- prediction(
  m3_anxiety, 
  at = list(anxiety_scale = c(0, 0.165, 0.33, 0.5, 0.66, 0.83, 1))
) %>% summary()

ggplot(m3_pred, aes(x = `at(anxiety_scale)`, y = Prediction, 
                    ymin = lower, ymax = upper)) +
  geom_point() +
  geom_ribbon(alpha = 0.4) +
  geom_line(alpha = 0.6, group = 1) +
  theme_custom() +
  labs(x = "COVID-related anxiety",
       y = "Predicted support for executive aggrandizement") +
  scale_y_continuous(limits = c(0.2, 0.6))

ggsave("figures/study 2/fig_2.png", height = 7, width = 9)
```

## Figure 3 

```{r}
# Manipulating data
exp1_lockdowns <- survey %>% 
  mutate(gov_group = strsplit(gov_group, "_") %>% sapply("[", 3)) %>% 
  dplyr::group_by(gov_group, covid_lockdowns) %>% 
  dplyr::summarise(mean_outcome = mean(exp1_outcome, na.rm = T),
                   lwr = lwr_conf(exp1_outcome),
                   upr = upr_conf(exp1_outcome)) %>% 
  mutate(label = paste0(round(mean_outcome,3), " [", round(lwr,3), ",", round(upr,3),"]"))

# results
ggplot(data = exp1_lockdowns,
       aes(x = covid_lockdowns, y = mean_outcome, 
           ymin = lwr, ymax = upr, col = gov_group, label = label,
           shape = gov_group, fill = gov_group)) +
  geom_point(size = 4, position = position_dodge(0.5)) +
  geom_errorbar(width = 0.15, size = 1, position = position_dodge(0.5)) +
  theme_custom() %+replace%
  theme(legend.position = "right") +
  scale_color_grey(name = "Premier's objective",
                     labels = c("Continue lockdown", "Open up")) +
  scale_fill_grey(name = "Premier's objective",
                     labels = c("Continue lockdown", "Open up")) +
  scale_shape_manual(values = c(21, 24),
                     name = "Premier's objective",
                     labels = c("Continue lockdown", "Open up")) +
  labs(y = "Proportion supporting shutdown",
       x = "Repondent preference on lockdowns") +
  scale_y_continuous(breaks = seq(0, 1, 0.1), limits = c(0, 1))

ggsave("figures/study 2/fig_3.png", height = 7, width = 11)
```

## Figure 4

```{r}
m4_pred <- prediction(
  m4_anxiety, 
  at = list(anxiety_scale = c(0, 0.165, 0.33, 0.5, 0.66, 0.83, 1),
            exp1_lockdown_congruence = c(0,1))
) %>% summary() %>% 
  mutate(`at(exp1_lockdown_congruence)` = factor(`at(exp1_lockdown_congruence)`))

ggplot(m4_pred, aes(x = `at(anxiety_scale)`, y = Prediction, 
                    ymin = lower, ymax = upper)) +
  geom_point(aes(col = `at(exp1_lockdown_congruence)`),
             size = 2) +
  geom_line(aes(linetype = `at(exp1_lockdown_congruence)`),
            alpha = 0.6) +
  geom_errorbar(aes(col = `at(exp1_lockdown_congruence)`),
                width = 0.01) +
  theme_custom(legend.position = "bottom") +
  scale_color_grey(end = 0.5,
                   name = "",
                   labels = c("Not congruent", "Congruent")) +
  scale_y_continuous(breaks = seq(0, 0.8, 0.2)) +
  scale_linetype_discrete(name = "",
                          labels = c("Not congruent", "Congruent")) +
  labs(x = "COVID-related anxiety",
       y = "Predicted agreement with shutting down")

ggsave("figures/study 2/fig_4.png", height = 7, width = 9)
```


# Study 3

## Time restrictions

```{r}
parliament.5s.task <- 
  filter(conjoint_parliament, task_time >= 5)

courts.5s.task <- 
  filter(conjoint_courts, task_time >= 5)
```

## Figure 5a

```{r}
cj_parliament <- prepare_conjoint(parliament.5s.task, "parliament")
plot_conjoint(cj_parliament, "partial") +
  scale_color_grey(end=0.6) +
  scale_x_continuous(limits = c(-.2, .3))

ggsave("figures/study 3/fig_5a.png", height = 6, width = 11)
```

## Figure 5b

```{r}
cj_courts <- prepare_conjoint(courts.5s.task, "courts")
plot_conjoint(cj_courts, "partial") +
  scale_color_grey(end=0.6) +
  scale_x_continuous(limits = c(-.2, .3))

ggsave("figures/study 3/fig_5b.png", height = 6, width = 11)
```

## Figure 6

```{r}
parl_anxiety <- 
  cregg::cj(
    filter(parliament.5s.task), 
    f1, by = ~anxiety_4cat, id = ~id
  ) %>% 
  filter(level == "Shut down Parliament") %>% 
  mutate(BY = factor(BY, levels = c("Never", "Rarely", "Sometimes", "Often")))

courts_anxiety <- 
  cregg::cj(
    filter(courts.5s.task), 
    f2, by = ~anxiety_4cat, id = ~id
  ) %>% 
  filter(level == "Ignore courts") %>% 
  mutate(BY = factor(BY, 
                              levels = c("Never", "Rarely", "Sometimes", "Often")))

anxiety_amce <- rbind(parl_anxiety, courts_anxiety)

ggplot(anxiety_amce,
       aes(x = estimate, y = fct_rev(level) , xmin = lower, xmax = upper,
           col = factor(BY), shape = factor(BY))) +
  geom_vline(xintercept = 0) +
  geom_point(size = 4, position = position_dodge(0.5)) +
  geom_errorbar(width = 0, size = 1, position = position_dodge(0.5),
                alpha = 0.8) +
  scale_color_grey(name = "Has felt anxiety...", end = 0.7) +
  scale_shape_manual(values = c(15, 16, 17, 18),
                     name = "Has felt anxiety...") +
  scale_x_continuous(breaks = c(-.2, -.1, 0, .1, .2, .3),
                           limits = c(-.20, .30)) +
  labs(x = "Estimated AMCE",
       y = "") +
  theme_cjoint() %+replace%
  theme(legend.position = "bottom",
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14))

ggsave("figures/study 3/fig_6.png", 
       height = 7, width = 12, dpi= 320)
```

### Formal test of AMCE differences
#### Parliament

```{r}
# Running regression by hand to test for significant difference in AMCEs by anxiety for the undemocratic attribute 
parl_anxious_reg <- 
  lm(selected ~ Gender*anxiety_4cat + Age*anxiety_4cat + 
       `Political experience`*anxiety_4cat +
       Party*anxiety_4cat + `Economic aid policy`*anxiety_4cat + 
       `Lockdown policy`*anxiety_4cat + `Legislative checks`*anxiety_4cat, 
     data = parliament.5s.task)

# Using clustered and robust SEs
parl_anxious_reg <- lmtest::coeftest(parl_anxious_reg,
                                     vcov = sandwich::vcovCL, cluster = ~id) 

# Extracting the coefficient and p-value for the interaction term:
# 'Often anxious':`Legislative checks`
# It represents the difference in AMCEs for the "Shut down legislature"
# attribute level beween the reference level of "anxiety_4cat" (Never anxious)
# and the often anxious
# p-value shows a significant different in AMCEs
parl_anxious_reg[rownames(parl_anxious_reg)=="anxiety_4catOften:`Legislative checks`Shut down Parliament",]
```

#### Courts

```{r}
# Running regression by hand to test for significant difference in AMCEs by anxiety for the undemocratic attribute 
courts_anxious_reg <- 
  lm(selected ~ Gender*anxiety_4cat + Age*anxiety_4cat + 
       `Political experience`*anxiety_4cat +
       Party*anxiety_4cat + `Economic aid policy`*anxiety_4cat + 
       `Lockdown policy`*anxiety_4cat + `Judicial checks`*anxiety_4cat, 
     data = courts.5s.task)

# Using clustered and robust SEs
courts_anxious_reg <- lmtest::coeftest(courts_anxious_reg,
                                       vcov = sandwich::vcovCL, cluster = ~id) 

# Extracting the coefficient and p-value for the interaction term:
# 'Often anxious':`Judicial checks`
# It represents the difference in AMCEs for the "Ignore courts"
# attribute level beween the reference level of "anxiety_4cat" (Never anxious)
# and the often anxious
# p-value shows no significant different in AMCEs
courts_anxious_reg[rownames(courts_anxious_reg)=="anxiety_4catOften:`Judicial checks`Ignore courts",]
```


## Figure 7a

```{r}
# Excluding tasks where one candidate had "lockdown until vaccine"
# and other one had "lockdown until fewer deaths" 
# i.e. restricting to tasks with stark choice on lockdowns 
conjoint_parliament_restrict <- parliament.5s.task %>% 
  dplyr::select(id:`Legislative checks`, covid_lockdowns, partyid,
                anxiety_4cat) %>% 
  pivot_wider(id_cols = c(id, task, covid_lockdowns, profile,
                          partyid, anxiety_4cat),
              names_sep = "_",
              names_from = profile,
              values_from = selected:`Legislative checks`) %>% 
  filter(!(`Lockdown policy_1` == "Lockdown until vaccine" & 
             `Lockdown policy_2` == "Lockdown until fewer deaths") &
           !(`Lockdown policy_1` == "Lockdown until fewer deaths" & 
             `Lockdown policy_2` == "Lockdown until vaccine")) %>% 
  pivot_longer(cols = selected_1:`Legislative checks_2`,
               names_to = c(".value", "profile"),
               names_pattern = "(.+)_(.+)")

# % voting for D- candidates under different scenarios of 
# congruence/divergence on lockdowns
defection_lockdown_parliament <- 
  conjoint_parliament_restrict %>% 
  conjoint_scenario(condition = "parliament", defect_cause = "lockdown")

# Plotting
defection_lockdown_parliament %>% 
  ggplot(aes(x = group, y = mean, ymin = lwr, ymax = upr, 
             label = label)) +
  geom_point(position = position_dodge(0.5), size = 3) +
  geom_errorbar(position = position_dodge(0.5), width = 0.1, size = 1,
                alpha = 0.8) +
  theme_custom() +
  labs(x = "Candidate-respondent congruence on lockdown policy",
       y = "Proportion choosing the undemocratic candidate") +
  geom_hline(yintercept = 0.5, linetype = "dashed") + 
  guides(col = FALSE) +
  scale_y_continuous(breaks = seq(0.1, 0.8, 0.1),
                     limits = c(0.11, 0.81)) +
  scale_x_discrete(
    labels = paste(defection_lockdown_parliament$group,
                   paste("n=", defection_lockdown_parliament$n, sep=""),
                   sep = "\n")
    )
ggsave("figures/study 3/fig_7a.png", height = 7, width = 10)
```

## Figure 7b 

```{r}
# Excluding tasks where one candidate had "lockdown until vaccine"
# and other one had "lockdown until fewer deaths" 
# i.e. restricting to tasks with stark choice on lockdowns 
conjoint_courts_restrict <- courts.5s.task %>%  
  dplyr::select(id:`Judicial checks`, covid_lockdowns, partyid,
                anxiety_4cat) %>% 
  pivot_wider(id_cols = c(id, task, covid_lockdowns, profile,
                          partyid, anxiety_4cat),
              names_sep = "_",
              names_from = profile,
              values_from = selected:`Judicial checks`) %>% 
  filter(!(`Lockdown policy_1` == "Lockdown until vaccine" & 
             `Lockdown policy_2` == "Lockdown until fewer deaths") &
           !(`Lockdown policy_1` == "Lockdown until fewer deaths" & 
             `Lockdown policy_2` == "Lockdown until vaccine")) %>% 
  pivot_longer(cols = selected_1:`Judicial checks_2`,
               names_to = c(".value", "profile"),
               names_pattern = "(.+)_(.+)")
  
# % voting for D- candidates under different scenarios of 
# congruence/divergence on lockdowns
defection_lockdown_courts <- 
  conjoint_courts_restrict %>% 
  conjoint_scenario(condition = "courts", defect_cause = "lockdown")

# Plotting
defection_lockdown_courts %>% 
  ggplot(aes(x = group, y = mean, ymin = lwr, ymax = upr, 
             label = label)) +
  geom_point(position = position_dodge(0.5), size = 3) +
  geom_errorbar(position = position_dodge(0.5), width = 0.1, size = 1,
                alpha = 0.8) +
  theme_custom() +
  labs(x = "Candidate-respondent congruence on lockdown policy",
       y = "Proportion choosing the undemocratic candidate") +
  geom_hline(yintercept = 0.5, linetype = "dashed") + 
  scale_y_continuous(breaks = seq(0.1, 0.8, 0.1),
                     limits = c(0.11, 0.81)) +
  scale_x_discrete(
    labels = paste(defection_lockdown_courts$group,
                   paste("n=", defection_lockdown_courts$n, sep=""),
                   sep = "\n")
    )

ggsave("figures/study 3/fig_7b.png", height = 7, width = 10)
```

## Figure 8a

```{r}
filter(conjoint_parliament_restrict, !is.na(anxiety_4cat)) %>% 
  conjoint_scenario(condition = "parliament",
                    defect_cause = "lockdown",
                    by = "anxiety_4cat") %>% 
  mutate(by = factor(by, levels = c("Never", "Rarely", "Sometimes", "Often"))) %>% 
  ggplot(aes(x = by, y = mean, ymin = lwr, ymax = upr, col = group,
             shape = group, fill = group)) +
  geom_point(position = position_dodge(0.5),
             size = 2.5) +
  geom_errorbar(position = position_dodge(0.5), width = 0) +
  scale_color_grey(name = "The respondent agrees with the lockdown stance of...",
                   start = 0.1, end = 0.7) +
  scale_fill_grey(name = "The respondent agrees with the lockdown stance of...",
                  start = 0.1, end = 0.7) +
  scale_shape_manual(values = 21:25,
                     name = "The respondent agrees with the lockdown stance of...") +
  labs(x = "Covid-related anxiety",
       name = "The respondent agrees with the lockdown stance of...",
       y = "Proportion choosing the undemocratic candidate") +
  theme_custom(base_size = 13) %+replace%
  theme(legend.position = "right")

ggsave("figures/study 3/fig_8a.png", height = 7, width = 10)
```

## Figure 8b

```{r}
filter(conjoint_courts_restrict, !is.na(anxiety_4cat)) %>% 
  conjoint_scenario(condition = "courts",
                    defect_cause = "lockdown",
                    by = "anxiety_4cat") %>% 
  mutate(by = factor(by, levels = c("Never", "Rarely", "Sometimes", "Often"))) %>% 
  ggplot(aes(x = by, y = mean, ymin = lwr, ymax = upr, col = group,
             shape = group, fill = group)) +
  geom_point(position = position_dodge(0.5),
             size = 2.5) +
  geom_errorbar(position = position_dodge(0.5), width = 0) +
  scale_color_grey(name = "The respondent agrees with the lockdown stance of...",
                   start = 0.1, end = 0.7) +
  scale_fill_grey(name = "The respondent agrees with the lockdown stance of...",
                  start = 0.1, end = 0.7) +
  scale_shape_manual(values = 21:25,
                     name = "The respondent agrees with the lockdown stance of...") +
  labs(x = "Covid-related anxiety",
       name = "The respondent agrees with the lockdown stance of...",
       y = "Proportion choosing the undemocratic candidate") +
  theme_custom(base_size = 13) %+replace%
  theme(legend.position = "right")

ggsave("figures/study 3/fig_8b.png", height = 7, width = 10)
```

### Formal test of difference

```{r}
conjoint_courts_scenario <- 
  conjoint_courts_restrict %>% 
  group_by(task, id) %>% 
  filter(length(unique(`Judicial checks`)) == 2) %>% 
  mutate(congruence = ifelse(covid_lockdowns == `Lockdown policy`,1,0),
         congruence_length = length(unique(congruence))) %>% 
  filter(`Judicial checks` == "Ignore courts") %>% 
  mutate(
    group = case_when(
      congruence == FALSE & congruence_length == 1 ~ "Neither candidate",
      congruence == TRUE & congruence_length == 2 ~ "Only the undemocratic candidate",
      congruence == FALSE & congruence_length == 2 ~ "Only the democratic candidate",
      congruence == TRUE & congruence_length == 1 ~ "Both candidates"
    )
  )
  
t.test(
  conjoint_courts_scenario$selected[conjoint_courts_scenario$anxiety_4cat == "Often" & conjoint_courts_scenario$group == "Only the undemocratic candidate"],
  conjoint_courts_scenario$selected[conjoint_courts_scenario$anxiety_4cat == "Never" & conjoint_courts_scenario$group == "Only the undemocratic candidate"]
)
```

# Appendix

## Figure A2a

```{r}
# Formal test: p = 0.39
cj_anova(data = parliament.5s.task,
         formula = f1,
         id = ~id,
         by = ~task)

# Plot of task-specific AMCEs 
parliament_tasks <- 
  cregg::cj(data = parliament.5s.task,
            formula = f1,
            id = ~id,
            level_order = "descending",
            by = ~task) %>% 
  mutate(estimate = ifelse(estimate == 0, NA, estimate))

ggplot(data = filter(parliament_tasks, feature != "Age"), 
       aes(x = estimate, y = level, col = task,
           xmin = lower, xmax = upper, shape = task)) +
  geom_vline(xintercept = 0) +
  geom_point(size = 2.5, position = position_dodge(0.8)) +
  geom_errorbar(width = 0, size = 1, position = position_dodge(0.8)) +
  labs(y = "") +
  theme_bw() +
  theme(legend.position = "right",
        text = element_text(size = 14)) +
  labs(title = "",
       x = "Estimated AMCE") +
  scale_color_grey(start = 0, end = 0.7) +
  scale_shape_manual(values = c(16, 15, 17, 10, 12, 8))

ggsave("figures/study 3/fig_a2a.png", height = 9, width = 8)
```

## Figure A2b

```{r}
# Formal test: p = 0.88
cj_anova(data = courts.5s.task,
         formula = f2,
         id = ~id,
         by = ~task)

# Plot of task-specific AMCEs 
courts_tasks <- 
  cregg::cj(data = courts.5s.task,
            formula = f2,
            id = ~id,
            level_order = "descending",
            by = ~task) %>% 
  mutate(estimate = ifelse(estimate == 0, NA, estimate))

ggplot(data = filter(courts_tasks, feature != "Age"), 
       aes(x = estimate, y = level, col = task, shape = task, 
           xmin = lower, xmax = upper)) +
  geom_vline(xintercept = 0) +
  geom_point(size = 2.5, position = position_dodge(0.8)) +
  geom_errorbar(width = 0, size = 1, position = position_dodge(0.8)) +
  labs(y = "") +
  theme_bw() +
  theme(legend.position = "right",
        text = element_text(size = 14)) +
  labs(title = "",
       x = "Estimated AMCE") +
  scale_color_grey(start = 0, end = 0.7) +
  scale_shape_manual(values = c(16, 15, 17, 10, 12, 8))
  

ggsave("figures/study 3/fig_a2b.png", height = 9, width = 8)
```

## Figure A3a

```{r}
# Formal test: p = 0.55
cj_anova(data = parliament.5s.task, 
         formula = f1,
         id = ~id,
         by = ~profile)

# Plot of task-specific AMCEs 
parliament_profiles <- 
  cregg::cj(data = parliament.5s.task,
            formula = f1,
            id = ~id,
            estimate = "mm",
            level_order = "descending",
            by = ~profile)

ggplot(data = filter(parliament_profiles, feature != "Age"), 
       aes(x = estimate, y = level, col = profile, 
           xmin = lower, xmax = upper, shape = profile)) +
  geom_vline(xintercept = 0.5) +
  geom_point(size = 3, position = position_dodge(0.7)) +
  geom_errorbar(width = 0, size = 1, position = position_dodge(0.7)) +
  labs(y = "") +
  theme_bw() +
  theme(legend.position = "right",
        text = element_text(size = 18)) +
  labs(title = "",
       x = "Estimated marginal mean") +
  scale_color_grey()

ggsave("figures/study 3/fig_a3a.png", height = 9, width = 11)
```

## Figure A3b

```{r}
# Formal test: p = 0.38
cj_anova(data = courts.5s.task, 
         formula = f2,
         id = ~id,
         by = ~profile)

# Plot of task-specific AMCEs 
courts_profiles <- 
  cregg::cj(data = courts.5s.task,
            formula = f2,
            id = ~id,
            estimate = "mm",
            level_order = "descending",
            by = ~profile)
  
ggplot(data = filter(courts_profiles, feature != "Age"), 
       aes(x = estimate, y = level, col = profile, 
           xmin = lower, xmax = upper, shape = profile)) +
  geom_vline(xintercept = 0.5) +
  geom_point(size = 3, position = position_dodge(0.7)) +
  geom_errorbar(width = 0, size = 1, position = position_dodge(0.7)) +
  labs(y = "") +
  theme_bw() +
  theme(legend.position = "right",
        text = element_text(size = 18)) +
  labs(title = "",
       x = "Estimated marginal mean") +
  scale_color_grey()

ggsave("figures/study 3/fig_a3b.png", height = 9, width = 11)
```

## Figure A4a

```{r}
# Formal test: p = 0.46
cj_anova(data = parliament.5s.task,
         formula = f1,
         id = ~id,
         by = ~democracy_row)

# Plot of task-specific AMCEs 
parliament_row <- 
  cregg::cj(data = parliament.5s.task,
            formula = f1,
            id = ~id,
            level_order = "descending",
            by = ~democracy_row) %>% 
  mutate(estimate = ifelse(estimate == 0, NA, estimate))

ggplot(data = filter(parliament_row, feature != "Age"), 
       aes(x = estimate, y = level, col = democracy_row, 
           xmin = lower, xmax = upper, shape = democracy_row)) +
  geom_vline(xintercept = 0) +
  geom_point(size = 3, position = position_dodge(0.7)) +
  geom_errorbar(width = 0, size = 1, position = position_dodge(0.7)) +
  labs(y = "") +
  theme_bw() +
  theme(legend.position = "right",
        text = element_text(size = 18)) +
  labs(title = "",
       x = "Estimated AMCE") +
  scale_color_grey(name = "Democracy row #") +
  scale_shape_manual(name = "Democracy row #",
                     values = c(16, 15, 17))

ggsave("figures/study 3/fig_a4a.png",
       height = 9, width = 11)
```
## Figure A4b

```{r}
# Formal test: p = 0.013
cj_anova(data = courts.5s.task, 
         formula = f2,
         id = ~id,
         by = ~democracy_row)

# Plot of task-specific AMCEs 
courts_row <- 
  cregg::cj(data = courts.5s.task, 
            formula = f2,
            id = ~id,
            estimate = "amce",
            level_order = "descending",
            by = ~democracy_row) %>% 
  mutate(estimate = ifelse(estimate == 0, NA, estimate))

ggplot(data = filter(courts_row, feature != "Age"), 
       aes(x = estimate, y = level, col = democracy_row, 
           xmin = lower, xmax = upper, shape = democracy_row)) +
  geom_vline(xintercept = 0) +
  geom_point(size = 3, position = position_dodge(0.7)) +
  geom_errorbar(width = 0, size = 1, position = position_dodge(0.7)) +
  labs(y = "") +
  theme_bw() +
  theme(legend.position = "right",
        text = element_text(size = 18)) +
  labs(title = "",
       x = "Estimated AMCE") +
  scale_color_grey(name = "Democracy row #") +
  scale_shape_manual(name = "Democracy row #",
                     values = c(16, 15, 17))

ggsave("figures/study 3/fig_a4b.png", height = 9, width = 11)
```



# Supplemental appendix

## Table S1

```{r}
survey_bal <- survey %>% 
  mutate(age = 2020 - birth_year) %>% 
  dplyr::select(
    age, gender, educ_4cat, partyid, deprivation_scale, 
    anxiety_scale, exp1_condition
  ) %>% 
  dummy_cols(
    remove_selected_columns = TRUE,
    select_columns = c("gender", "educ_4cat", "partyid")
    ) %>% 
  dplyr::select(exp1_condition, age, deprivation_scale, gender_Male,
                `educ_4cat_College graduate`) %>% 
  mutate(exp1_condition = recode(exp1_condition,
                                 "New Democrat/lockdown" = "NDP/lockdown",
                                 "New Democrat/open" = "NDP/open")) %>% 
  dplyr::rename(`Condition` = exp1_condition,
                Age = age,
                `Deprivation` = deprivation_scale,
                Male = gender_Male,
                `College Graduate` = `educ_4cat_College graduate`) 

bal_table <- survey_bal %>%   
  group_by(`Condition`) %>% 
  summarise(across(everything(), mean, na.rm = TRUE)) %>% 
  mutate(across(where(is.numeric), round, 2),
         Condition = as.character(Condition))

# extracting f-stats: this makes a list of regression output
bal_regs <- lm(cbind(Age, Deprivation, Male, `College Graduate`) ~ Condition,
   data = survey_bal) %>% 
  summary()

# for each element in the list, the 10th sub-element is the fstat
fstats <- lapply(bro, "[", 10) %>% lapply("[[", 1) %>% lapply("[", 1) %>% 
  unlist() %>% round(3) %>% paste("F = ", ., sep = "")

# final table 
rbind(bal_table, c("", fstats))
```


## Figure S1

```{r}
# Manipulating data
exp1_overall <- filter(survey, !is.na(gov_pid)) %>% 
  mutate(gov_group = strsplit(gov_group, "_") %>% sapply("[", 3)) %>% 
  dplyr::group_by(gov_group, gov_pid) %>% 
  dplyr::summarise(mean_outcome = mean(exp1_outcome, na.rm = T),
                   lwr = lwr_conf(exp1_outcome),
                   upr = upr_conf(exp1_outcome)) %>% 
  mutate(label = paste0(round(mean_outcome,3), " [", round(lwr,3), ",", round(upr,3),"]"))
  

# Overall results
ggplot(data = exp1_overall,
       aes(x = gov_pid, y = mean_outcome, ymin = lwr, ymax = upr, 
           label = label, shape = gov_group, col = gov_group)) +
  geom_point(size = 3.5, position = position_dodge(0.5)) +
  geom_errorbar(width = 0.15, size = 1, position = position_dodge(0.5)) +
  theme_custom() %+replace%
  theme(legend.position = "right") +
  scale_color_grey(end = 0.7,
                   name = "Premier's objective",
                     labels = c("Continue lockdown", "Open up")) +
  scale_shape_manual(values = c(16, 15),
                     name = "Premier's objective",
                     labels = c("Continue lockdown", "Open up")) +
  labs(y = "Proportion supporting shutdown",
       x = "Partisan affiliation of the premier") +
  scale_y_continuous(breaks = seq(0, 1, 0.1), limits = c(0, 1))

ggsave("figures/study 2/fig_s1.png.png", height = 7, width = 10)
```

## Figure S2

```{r}
exp1_pid_congruence_pid <- 
  filter(survey, !is.na(gov_pid) & 
           partyid %in% c("Conservative", "Liberal", "NDP")) %>% 
  mutate(gov_group = strsplit(gov_group, "_") %>% sapply("[", 3)) %>% 
  dplyr::group_by(partyid, gov_pid) %>% 
  dplyr::summarise(mean_outcome = mean(exp1_outcome, na.rm = T),
                   lwr = lwr_conf(exp1_outcome),
                   upr = upr_conf(exp1_outcome)) 

# Plotting
ggplot(data = exp1_pid_congruence_pid,
       aes(x = partyid, y = mean_outcome, ymin = lwr, ymax = upr, 
           col = gov_pid, shape = gov_pid)) +
  geom_point(size = 3.5, position = position_dodge(0.5)) +
  geom_errorbar(width = 0.15, size = 1, position = position_dodge(0.5)) +
  theme_custom()  %+replace%
  theme(legend.position = "right") +
  labs(y = "Proportion supporting shutdown",
       x = "Respondent's party ID") +
  scale_color_manual(values =  c("#1A4782","#D71920", "#F37021", "grey"),
                     name = "Premier's party",
                     labels = c("Conservative", "Liberal", 
                                "New Democrat", "No party")) +
  scale_shape_manual(values = c(16, 15, 17, 18),
                     name = "Premier's party",
                     labels = c("Conservative", "Liberal", 
                                "New Democrat", "No party")) +
  scale_y_continuous(breaks = seq(0, 1, 0.1), limits = c(0, 1))

ggsave("figures/study 2/fig_s2.png", 
       height = 7, width = 10)

# T-test: Conservative respondents voting for C and L candidate
t.test(survey$exp1_outcome[survey$partyid=="Conservative" &
                             survey$gov_pid=="Conservative"],
       survey$exp1_outcome[survey$partyid=="Conservative" &
                             survey$gov_pid=="Liberal"])
# T-test: NDP respondents voting for C and L candidate
t.test(survey$exp1_outcome[survey$partyid=="NDP" &
                             survey$gov_pid=="Conservative"],
       survey$exp1_outcome[survey$partyid=="NDP" &
                             survey$gov_pid=="Liberal"])
```

## Figure S3

```{r}
dem_estimate <- data.frame(restriction = seq(0, 30, .25),
                           estimate_parliament = NA,
                           lwr_parliament = NA,
                           upr_parliament = NA,
                           estimate_courts = NA,
                           lwr_courts = NA,
                           upr_courts = NA)

for(i in 0:(nrow(dem_estimate)-1)){
  conjoint_restricted_parliament <- filter(conjoint_parliament, task_time > i*0.25)
  dem_estimate[i+1,2:4] <- cregg::cj(conjoint_restricted_parliament,
            formula = f1,
            id = ~id) %>% 
    filter(level == "Shut down Parliament") %>% 
    dplyr::select(estimate, lower, upper) %>% 
    as.vector() %>% 
    as.numeric()
  
  conjoint_restricted_courts <- filter(conjoint_courts, task_time > i*0.25)
  dem_estimate[i+1,5:7] <- cregg::cj(conjoint_restricted_courts,
            formula = f2,
            id = ~id) %>% 
    filter(level == "Ignore courts") %>% 
    dplyr::select(estimate, lower, upper) %>% 
    as.vector() %>% 
    as.numeric()
}

restrict_parliament <- ggplot(data = dem_estimate,
       aes(x = restriction, y = estimate_parliament, ymin = lwr_parliament, ymax = upr_parliament)) +
  geom_line(aes(x = restriction, y = estimate_parliament)) +
  geom_ribbon(aes(x = restriction, y = estimate_parliament, ymin = lwr_parliament, ymax = upr_parliament), 
              alpha = 0.4) +
  labs(y = "AMCE for shutting down parliament",
       x = "Time restriction (seconds)") +
  theme_custom() +
scale_y_continuous(limits = c(-0.12, -0.03),
                   breaks = seq(-0.12, -0.03, 0.01))

restrict_courts <- ggplot(data = dem_estimate) +
  geom_line(aes(x = restriction, y = estimate_courts)) +
  geom_ribbon(aes(x = restriction, y = estimate_courts, ymin = lwr_courts, ymax = upr_courts), 
              alpha = 0.4) +
  labs(y = "AMCE for ignoring courts",
       x = "Time restriction (seconds)") +
  theme_custom()+
scale_y_continuous(limits = c(-0.12, -0.03),
                   breaks = seq(-0.12, -0.03, 0.01))

ggarrange(restrict_parliament, restrict_courts)
ggsave("figures/study 3/fig_s3.png", height = 7, width = 12)
```

## Figure S4a

```{r}
cj_parliament <- prepare_conjoint(parliament.5s.task, "parliament")
plot_conjoint(cj_parliament)
ggsave("figures/study 3/fig_s4a.png", height = 9, width = 9)
```

## Figure S4b

```{r}
cj_courts <- prepare_conjoint(courts.5s.task, "courts")
plot_conjoint(cj_courts)
ggsave("figures/study 3/fig_s4b.png", height = 9, width = 9)
```

## Figure S5a

```{r}
parl_anxious1 <- 
  cregg::cj(
    parliament.5s.task, 
    f1, by = ~anxious1_think, id = ~id
  ) %>% 
  filter(level == "Shut down Parliament") %>% 
  mutate(BY = as.character(BY) %>% 
           dplyr::recode("0" = "Never", "0.33" = "Rarely", 
                         "0.66" = "Sometimes", "1" = "Often") %>% 
           factor(levels = c("Never", "Rarely", "Sometimes", "Often")))

courts_anxious1 <- 
  cregg::cj(
    courts.5s.task, 
    f2, by = ~anxious1_think, id = ~id
  ) %>% 
  filter(level == "Ignore courts") %>% 
  mutate(BY = as.character(BY) %>% 
           dplyr::recode("0" = "Never", "0.33" = "Rarely", 
                         "0.66" = "Sometimes", "1" = "Often") %>% 
           factor(levels = c("Never", "Rarely", "Sometimes", "Often")))

anxious1_amce <- rbind(parl_anxious1, courts_anxious1)

ggplot(anxious1_amce,
       aes(x = estimate, y = fct_rev(level) , xmin = lower, xmax = upper,
           col = factor(BY), shape = factor(BY))) +
  geom_vline(xintercept = 0) +
  geom_point(size = 4, position = position_dodge(0.5)) +
  geom_errorbar(width = 0, size = 1, position = position_dodge(0.5),
                alpha = 0.8) +
  scale_color_grey(name = "Has felt anxiety...", end = 0.7) +
  scale_shape_manual(values = c(15, 16, 17, 18),
                     name = "Has felt anxiety...") +
  scale_x_continuous(breaks = c(-.2, -.1, 0, .1, .2, .3),
                           limits = c(-.2, .3)) +
  labs(x = "Estimated AMCE",
       y = "") +
  theme_cjoint() %+replace%
  theme(legend.position = "bottom",
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14))

ggsave("figures/study 3/fig_s5a.png", height = 7, width = 12)
```

## Figure S5b

```{r}
parl_anxious2 <- 
  cregg::cj(
    parliament.5s.task, 
    f1, by = ~anxious2_activities, id = ~id
  ) %>% 
  filter(level == "Shut down Parliament") %>% 
  mutate(BY = as.character(BY) %>% 
           dplyr::recode("0" = "Never", "0.33" = "Rarely", 
                         "0.66" = "Sometimes", "1" = "Often") %>% 
           factor(levels = c("Never", "Rarely", "Sometimes", "Often")))

courts_anxious2 <- 
  cregg::cj(
    courts.5s.task, 
    f2, by = ~anxious2_activities, id = ~id
  ) %>% 
  filter(level == "Ignore courts") %>% 
  mutate(BY = as.character(BY) %>% 
           dplyr::recode("0" = "Never", "0.33" = "Rarely", 
                         "0.66" = "Sometimes", "1" = "Often") %>% 
           factor(levels = c("Never", "Rarely", "Sometimes", "Often")))

anxious2_amce <- rbind(parl_anxious2, courts_anxious2)

ggplot(anxious2_amce,
       aes(x = estimate, y = fct_rev(level) , xmin = lower, xmax = upper,
           col = factor(BY), shape = factor(BY))) +
  geom_vline(xintercept = 0) +
  geom_point(size = 4, position = position_dodge(0.5)) +
  geom_errorbar(width = 0, size = 1, position = position_dodge(0.5),
                alpha = 0.8) +
  scale_color_grey(name = "Has felt anxiety...", end = 0.7) +
  scale_shape_manual(values = c(15, 16, 17, 18),
                     name = "Has felt anxiety...") +
  scale_x_continuous(breaks = c(-.2, -.1, 0, .1, .2, .3),
                           limits = c(-.2, .3)) +
  labs(x = "Estimated AMCE",
       y = "") +
  theme_cjoint() %+replace%
  theme(legend.position = "bottom",
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14))

ggsave("figures/study 3/fig_s5b.png", height = 7, width = 12)
```







